Luis Alejandro RodrÃguez Arenas Cod. 202321287¶
BLOQUE 1 Imports y configuración¶
In [1]:
import numpy as np
from numpy import linspace
import matplotlib.pyplot as plt
from sympy import symbols, sin, cos, pi, integrate, lambdify, simplify, latex, exp
from sympy import init_printing
from sympy import Piecewise
from matplotlib.animation import FuncAnimation
init_printing(use_latex='mathjax')
In [2]:
x, t = symbols('x t', real=True)
BLOQUE 2 - Parámetros fÃsicos del problema¶
ej 1 10-5¶
In [ ]:
L = 50
alpha = 2
T0 = 0
TL = 0
f = 20
print("Función inicial f(x) =")
display(f)
ej 1 10-6¶
In [3]:
L = 30
alpha = 2
T0 = 20
TL = 50
f = 60-2*x
print("Función inicial f(x) =")
display(f)
Función inicial f(x) =
$\displaystyle 60 - 2 x$
ej 1 10-7¶
In [16]:
L = 30
a = 4
T0 = 0
TL = 0
f = Piecewise(
(x/10, (0 <= x) & (x <= 10)),
((30-x)/20, (10 < x) & (x <= 30)),
(0, x > 30),
)
print("Función inicial f(x) por tramos =")
display(f)
Función inicial f(x) por tramos =
$\displaystyle \begin{cases} \frac{x}{10} & \text{for}\: x \geq 0 \wedge x \leq 10 \\\frac{3}{2} - \frac{x}{20} & \text{for}\: x \leq 30 \wedge x > 10 \\0 & \text{for}\: x > 30 \end{cases}$
In [ ]:
L = 30
alpha = 2
T0 = 0
TL = 0
f = sin(15*pi*x/L) + sin(35*pi*x/L)
Separación en parte estacionaria + transitoria¶
In [17]:
# Parte estacionaria v(x): resolviendo v''(x)=0, v(0)=T0, v(L)=TL
v = T0 + (TL - T0) * (x/L)
print("Solución estacionaria v(x) =")
display(v)
# Condición inicial para v(x,0): v0(x) = f(x) – s(x)
v0 = simplify(f - v)
print("\nv(x,0) = f(x) - s(x) =")
display(v0)
Solución estacionaria v(x) =
$\displaystyle 0$
v(x,0) = f(x) - s(x) =
$\displaystyle \begin{cases} \frac{x}{10} & \text{for}\: x \geq 0 \wedge x \leq 10 \\\frac{3}{2} - \frac{x}{20} & \text{for}\: x \leq 30 \wedge x > 10 \\0 & \text{for}\: x > 30 \\\text{NaN} & \text{otherwise} \end{cases}$
BLOQUE 4 - Serie de Fourier¶
In [18]:
# Fórmula simbólica del coeficiente c_n
n = symbols('n', integer=True, positive=True)
cn_formula = (2/L) * integrate(v0 * sin(n*pi*x/L), (x, 0, L))
print("Coeficiente general c_n =")
display(simplify(cn_formula))
Coeficiente general c_n =
$\displaystyle \frac{9.0 \sin{\left(\frac{\pi n}{3} \right)}}{\pi^{2} n^{2}}$
In [6]:
# Para calor
print("u(x,T) =")
display(cn_formula * exp(-(alpha*n*pi/L)**2 * t) * sin(n*pi*x/L))
u(x,T) =
$\displaystyle \left(\frac{100.0 \left(-1\right)^{n}}{\pi n} + \frac{80.0}{\pi n}\right) e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}$
In [19]:
# Para onda
print("u(x,t) =")
display(cn_formula * cos(-(alpha*n*pi/L) * t) * sin(n*pi*x/L))
u(x,t) =
$\displaystyle \frac{9.0 \sin{\left(\frac{\pi n}{3} \right)} \sin{\left(\frac{\pi n x}{30} \right)} \cos{\left(\frac{\pi n t}{15} \right)}}{\pi^{2} n^{2}}$
Cálculo de coeficientes numéricos (Calor)¶
In [7]:
N = 100
cn = []
for k in range(1, N+1):
val = float(cn_formula.subs(n, k))
cn.append(val)
print("\nPrimeros coeficientes c_n (numéricos):")
for i, c in enumerate(cn[:8], 1):
print(f"c{i} = {c:.6f}")
u_series_sym = 0
for k in range(1, N+1):
ck = cn_formula.subs(n, k)
term = ck * exp(-(alpha*n*pi/L)**2 * t) * sin(n*pi*x/L)
u_series_sym += term
print("\nSolución v(x,t) truncada simbólica (antes de sumar v(x)):")
display(simplify(u_series_sym))
u_full_sym = v + u_series_sym
print("\nSolución completa u(x,t)")
display(simplify(u_full_sym))
Primeros coeficientes c_n (numéricos): c1 = -6.366198 c2 = 28.647890 c3 = -2.122066 c4 = 14.323945 c5 = -1.273240 c6 = 9.549297 c7 = -0.909457 c8 = 7.161972 Solución v(x,t) truncada simbólica (antes de sumar v(x)):
$\displaystyle \frac{346.17298348015 e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}}{\pi}$
Solución completa u(x,t)
$\displaystyle x + 20 + \frac{346.17298348015 e^{- \frac{\pi^{2} n^{2} t}{225}} \sin{\left(\frac{\pi n x}{30} \right)}}{\pi}$
Cálculo de coeficientes numéricos (Onda)¶
In [20]:
# Coeficiente correcto para ONDA
A_n_formula = (2/L) * integrate(f * sin(n*pi*x/L), (x, 0, L))
N = 100
A = []
for k in range(1, N+1):
A.append(float(A_n_formula.subs(n, k)))
print("Coeficientes A_n de la onda:")
for i, c in enumerate(A[:8], 1):
print(f"A{i} = {c:.6f}")
def u_series_onda(x_array, t_value, A, N):
total = np.zeros_like(x_array)
for k in range(1, N+1):
total += A[k-1] * np.sin(k*np.pi*x_array/L) * np.cos((k*np.pi*a/L)*t_value)
return total
def u_sol_onda(x_array, t_value):
return u_series_onda(x_array, t_value, A, N)
Coeficientes A_n de la onda: A1 = 0.789720 A2 = 0.197430 A3 = 0.000000 A4 = -0.049358 A5 = -0.031589 A6 = 0.000000 A7 = 0.016117 A8 = 0.012339
BLOQUE 5 - Construcción de u(x,t)¶
Calor¶
In [8]:
x_vals = np.linspace(0, L, 300)
def v_series_calor(x_array, t_value, cn, N):
total = np.zeros_like(x_array)
for k in range(1, N+1):
total += cn[k-1] * np.sin(k*np.pi*x_array/L) * np.exp(-alpha*(k*np.pi/L)**2 * t_value)
return total
v_func = lambdify(x, v, 'numpy')
def u_sol_calor(x_array, t_value):
return v_func(x_array) + v_series_calor(x_array, t_value, cn, N)
u_sol_calor(x_vals, 0)
Out[8]:
array([20. , 45.21579184, 61.84221647, 66.55056011, 62.68915318,
57.07310592, 54.92216581, 56.86623359, 59.89584709, 60.82157044,
59.04915746, 56.5541714 , 55.61603767, 56.6844463 , 58.25590197,
58.56542695, 57.28372131, 55.6454408 , 55.0841114 , 55.83566643,
56.83905447, 56.89287802, 55.83423964, 54.60252645, 54.22685952,
54.81180989, 55.50833404, 55.42075765, 54.49207593, 53.50127626,
53.23719968, 53.71732452, 54.21885641, 54.04075526, 53.19892705,
52.36984744, 52.18083512, 52.58739073, 52.95208813, 52.71062262,
51.93227901, 51.22093087, 51.08620101, 51.43709715, 51.69904853,
51.41050338, 50.68163564, 50.06106032, 49.96759107, 50.27397334,
50.45487023, 50.12984184, 49.44145415, 48.89393664, 48.83297236,
49.10219033, 49.21668866, 48.86251043, 48.20853564, 47.72181226,
47.68713577, 47.92424236, 47.9827078 , 47.60470862, 46.98090689,
46.54613737, 46.53313528, 46.7417106 , 46.75174329, 46.35395393,
45.75728636, 45.36788925, 45.37300994, 45.55564437, 45.52298103,
45.10855567, 44.53680701, 44.18775212, 44.20817365, 44.36676601,
44.29584122, 43.86732205, 43.31886259, 43.00622109, 43.03963833,
43.17558768, 43.06989791, 42.62938806, 42.10301784, 41.82366521,
41.86814803, 41.98248127, 41.84482917, 41.39410991, 40.88895354,
40.64036724, 40.69426294, 40.78772192, 40.62038499, 40.1609976 ,
39.67643169, 39.45654966, 39.5184138 , 39.59151611, 39.39636605,
38.92967062, 38.46527275, 38.27239229, 38.34093835, 38.39402021,
38.172609 , 37.69982784, 37.25534039, 37.08804454, 37.1621063 ,
37.19535311, 36.94897619, 36.47122657, 36.046531 , 35.90363422,
35.98213698, 35.99560492, 35.725348 , 35.24366769, 34.8387665 ,
34.71927417, 34.80121207, 34.79484297, 34.5016171 , 34.01698465,
33.63198922, 33.53506735, 33.61948496, 33.59311609, 33.27768389,
32.79103541, 32.42615846, 32.35111096, 32.43708793, 32.3904575 ,
32.05345266, 31.5656961 , 31.22124803, 31.16750005, 31.25413769,
31.18688673, 30.82882821, 30.34085603, 30.01724476, 29.98433075,
30.07073989, 29.98241075, 29.60371268, 29.11641364, 28.81414773,
28.80170349, 28.88699291, 28.7770245 , 28.37800228, 27.8922729 ,
27.61196805, 27.6197263 , 27.70299124, 27.57071073, 27.15158377,
26.66834008, 26.41072929, 26.43851859, 26.51882863, 26.36343927,
25.92433049, 25.44452061, 25.21046855, 25.25821546, 25.3346013 ,
25.15516562, 24.69609759, 24.22071575, 24.01123828, 24.07897303,
24.15041139, 23.94582874, 23.46671602, 22.99681899, 22.81310894,
22.90097513, 22.96637089, 22.7353478 , 22.23598499, 21.77271169,
21.61617286, 21.72444203, 21.78260639, 21.52361756, 21.00366187,
20.54825777, 20.42054978, 20.54964203, 20.59926517, 20.31050192,
19.7694488 , 19.32329682, 19.2263945 , 19.37690719, 19.4165231 ,
19.09582465, 18.53297408, 18.09763481, 18.03390789, 18.20665526,
18.23459547, 17.87935626, 17.29376619, 16.87103144, 16.84335279,
17.03942085, 17.05375213, 16.66079485, 16.05121621, 15.64318192,
15.65507726, 15.87590085, 15.87433931, 15.43973752, 14.80452246,
14.41369056, 14.46954978, 14.71702246, 14.69681203, 14.21563679,
13.55260642, 13.1820308 , 13.28741371, 13.56404803, 13.5217837 ,
12.98773191, 12.29398071, 11.94748305, 12.10957474, 12.41874273,
12.35010478, 11.75493651, 11.02653404, 10.70903413, 10.93734691,
11.28365344, 11.1829927 , 10.51564645, 9.74716512, 9.46520749,
9.77270872, 10.1625959 , 10.02225659, 9.2673933 , 8.45112487,
8.21376146, 8.61877974, 9.06155727, 8.87070887, 8.00617681,
7.13075527, 6.95111824, 7.48077522, 7.99049502, 7.7329729 ,
6.7250672 , 5.77286293, 5.67119607, 6.36810831, 6.96727509,
6.61721362, 5.41094316, 4.35262857, 4.36277122, 5.29964075,
6.0274404 , 5.5392927 , 4.03566136, 2.81724289, 3.00267863,
4.31934875, 5.25311901, 4.53440507, 2.52634234, 1.03131863,
1.5348144 , 3.55779259, 4.88545709, 3.69735932, 0.62219005,
-1.48410087, -0.21203353, 3.6298674 , 6.05004443, 3.34766987,
-3.62505196, -8.34057928, -2.34202258, 18.49211875, 50. ])
Onda¶
In [21]:
x_vals = np.linspace(0, L, 300)
def u_series_onda(x_array, t_value, A, N):
total = np.zeros_like(x_array)
for k in range(1, N+1):
total += A[k-1] * np.sin(k*np.pi*x_array/L) * np.cos((k*np.pi*a/L)*t_value)
return total
def u_sol_onda(x_array, t_value):
return u_series_onda(x_array, t_value, A, N)
u_sol_onda(x_vals, 0)
Out[21]:
array([0.00000000e+00, 1.00343941e-02, 2.00690687e-02, 3.01027554e-02,
4.01339423e-02, 5.01636934e-02, 6.01958188e-02, 7.02329527e-02,
8.02725683e-02, 9.03082994e-02, 1.00336495e-01, 1.10361270e-01,
1.20391570e-01, 1.30431962e-01, 1.40476607e-01, 1.50513945e-01,
1.60538600e-01, 1.70558272e-01, 1.80587183e-01, 1.90631421e-01,
2.00681262e-01, 2.10719770e-01, 2.20740249e-01, 2.30754609e-01,
2.40782575e-01, 2.50831347e-01, 2.60886641e-01, 2.70925863e-01,
2.80941413e-01, 2.90950150e-01, 3.00977646e-01, 3.11031783e-01,
3.21092900e-01, 3.31132335e-01, 3.41142033e-01, 3.51144712e-01,
3.61172268e-01, 3.71232804e-01, 3.81300260e-01, 3.91339335e-01,
4.01342015e-01, 4.11338027e-01, 4.21366273e-01, 4.31434532e-01,
4.41509041e-01, 4.51547061e-01, 4.61541205e-01, 4.71529706e-01,
4.81559421e-01, 4.91637162e-01, 5.01719725e-01, 5.11755801e-01,
5.21739354e-01, 5.31719152e-01, 5.41751362e-01, 5.51841011e-01,
5.61933058e-01, 5.71965986e-01, 5.81936050e-01, 5.91905415e-01,
6.01941558e-01, 6.12046621e-01, 6.22150267e-01, 6.32178300e-01,
6.42130570e-01, 6.52086887e-01, 6.62129129e-01, 6.72254978e-01,
6.82373501e-01, 6.92393894e-01, 7.02321538e-01, 7.12260628e-01,
7.22312536e-01, 7.32468056e-01, 7.42606892e-01, 7.52614871e-01,
7.62506024e-01, 7.72420643e-01, 7.82488836e-01, 7.92690392e-01,
8.02859483e-01, 8.12845460e-01, 8.22676624e-01, 8.32552620e-01,
8.42651779e-01, 8.52934910e-01, 8.63155356e-01, 8.73095134e-01,
8.82808667e-01, 8.92612615e-01, 9.02787010e-01, 9.13255031e-01,
9.23584207e-01, 9.33382600e-01, 9.42761832e-01, 9.52384821e-01,
9.62908697e-01, 9.74172342e-01, 9.84799512e-01, 9.92665184e-01,
9.96028152e-01, 9.94574765e-01, 9.89622574e-01, 9.83327655e-01,
9.77464303e-01, 9.72626103e-01, 9.68324356e-01, 9.63758414e-01,
9.58575439e-01, 9.53061104e-01, 9.47744872e-01, 9.42875010e-01,
9.38241973e-01, 9.33445966e-01, 9.28297652e-01, 9.22967190e-01,
9.17777975e-01, 9.12881337e-01, 9.08131037e-01, 9.03254207e-01,
8.98126194e-01, 8.92877234e-01, 8.87744197e-01, 8.82831893e-01,
8.78021509e-01, 8.73103193e-01, 8.67987539e-01, 8.62784821e-01,
8.57683811e-01, 8.52761552e-01, 8.47914387e-01, 8.42970492e-01,
8.37863006e-01, 8.32690319e-01, 8.27610293e-01, 8.22681225e-01,
8.17809012e-01, 8.12847523e-01, 8.07745794e-01, 8.02594371e-01,
7.97529366e-01, 7.92595370e-01, 7.87704831e-01, 7.82730322e-01,
7.77632859e-01, 7.72497435e-01, 7.67443863e-01, 7.62506136e-01,
7.57601473e-01, 7.52616811e-01, 7.47522639e-01, 7.42399810e-01,
7.37355343e-01, 7.32414686e-01, 7.27498694e-01, 7.22505797e-01,
7.17414254e-01, 7.12301697e-01, 7.07264732e-01, 7.02321701e-01,
6.97396329e-01, 6.92396549e-01, 6.87307168e-01, 6.82203230e-01,
6.77172618e-01, 6.72227605e-01, 6.67294265e-01, 6.62288592e-01,
6.57201039e-01, 6.52104505e-01, 6.47079388e-01, 6.42132678e-01,
6.37192424e-01, 6.32181605e-01, 6.27095636e-01, 6.22005589e-01,
6.16985310e-01, 6.12037112e-01, 6.07090747e-01, 6.02075363e-01,
5.96990800e-01, 5.91906532e-01, 5.86890573e-01, 5.81941041e-01,
5.76989191e-01, 5.71969706e-01, 5.66886416e-01, 5.61807371e-01,
5.56795317e-01, 5.51844563e-01, 5.46887724e-01, 5.41864513e-01,
5.36782400e-01, 5.31708134e-01, 5.26699644e-01, 5.21747751e-01,
5.16786321e-01, 5.11759695e-01, 5.06678689e-01, 5.01608844e-01,
4.96603633e-01, 4.91650660e-01, 4.86684963e-01, 4.81655182e-01,
4.76575233e-01, 4.71509517e-01, 4.66507345e-01, 4.61553333e-01,
4.56583634e-01, 4.51550920e-01, 4.46471996e-01, 4.41410167e-01,
4.36410830e-01, 4.31455804e-01, 4.26482322e-01, 4.21446866e-01,
4.16368948e-01, 4.11310806e-01, 4.06314124e-01, 4.01358099e-01,
3.96381017e-01, 3.91342985e-01, 3.86266064e-01, 3.81211442e-01,
3.76217260e-01, 3.71260239e-01, 3.66279710e-01, 3.61239250e-01,
3.56163325e-01, 3.51112084e-01, 3.46120263e-01, 3.41162242e-01,
3.36178394e-01, 3.31135637e-01, 3.26060716e-01, 3.21012737e-01,
3.16023153e-01, 3.11064123e-01, 3.06077065e-01, 3.01032127e-01,
2.95958223e-01, 2.90913407e-01, 2.85925949e-01, 2.80965894e-01,
2.75975715e-01, 2.70928704e-01, 2.65855835e-01, 2.60814098e-01,
2.55828665e-01, 2.50867563e-01, 2.45874342e-01, 2.40825354e-01,
2.35753543e-01, 2.30714816e-01, 2.25731314e-01, 2.20769141e-01,
2.15772942e-01, 2.10722065e-01, 2.05651339e-01, 2.00615563e-01,
1.95633908e-01, 1.90670635e-01, 1.85671509e-01, 1.80618826e-01,
1.75549217e-01, 1.70516343e-01, 1.65536456e-01, 1.60572049e-01,
1.55570043e-01, 1.50515629e-01, 1.45447171e-01, 1.40417159e-01,
1.35438967e-01, 1.30473391e-01, 1.25468538e-01, 1.20412465e-01,
1.15345194e-01, 1.10318015e-01, 1.05341448e-01, 1.00374665e-01,
9.53669936e-02, 9.03093269e-02, 8.52432843e-02, 8.02189119e-02,
7.52439069e-02, 7.02758747e-02, 6.52654056e-02, 6.02062077e-02,
5.51414362e-02, 5.01198537e-02, 4.51463503e-02, 4.01770244e-02,
3.51637715e-02, 3.01031008e-02, 2.50396463e-02, 2.00208429e-02,
1.50487844e-02, 1.00781174e-02, 5.06208853e-03, 2.67536267e-16])
BLOQUE 6 - Gráficas¶
Calor¶
In [9]:
x_vals = np.linspace(0, L, 300)
t_list = [0, 1, 10, 100, 500]
plt.figure(figsize=(8,5))
for ti in t_list:
ui = u_sol_calor(x_vals, ti)
plt.plot(x_vals, ui, label=f"t={ti}s")
plt.title("Solución analÃtica por Serie de Fourier")
plt.xlabel("x")
plt.ylabel("Temperatura u(x,t)")
plt.legend()
plt.grid(True)
plt.show()
Onda¶
In [ ]:
x_vals = np.linspace(0, L, 300)
t_list = [0, 1, 10, 100, 500]
plt.figure(figsize=(8,5))
for ti in t_list:
ui = u_sol_onda(x_vals, ti)
plt.plot(x_vals, ui, label=f"t={ti}s")
plt.title("Solución analÃtica por Serie de Fourier")
plt.xlabel("x")
plt.ylabel("Temperatura u(x,t)")
plt.legend()
plt.grid(True)
plt.show()
BLOQUE 7 - Animación¶
Calor¶
In [10]:
fig, ax = plt.subplots(figsize=(8,5))
line, = ax.plot([], [], lw=2)
ax.set_xlim(0, L)
ax.set_ylim(np.min(u_sol_calor(x_vals, 500)), np.max(u_sol_calor(x_vals, 0)))
ax.grid(True)
t_anim = np.linspace(0, 200, 300)
def init():
line.set_data([], [])
return line,
def update(frame):
t_val = t_anim[frame]
u_val = u_sol_calor(x_vals, t_val)
line.set_data(x_vals, u_val)
ax.set_title(f"t = {t_val:.1f} s")
return line,
anim = FuncAnimation(
fig,
update,
frames=len(t_anim),
init_func=init,
interval=100,
blit=False
)
from IPython.display import HTML
HTML(anim.to_jshtml())
Out[10]:
Onda¶
In [23]:
fig, ax = plt.subplots(figsize=(8,5))
line, = ax.plot([], [], lw=2)
modelo = "onda"
if modelo == "calor":
u_test_min = u_sol_calor(x_vals, 500)
u_test_max = u_sol_calor(x_vals, 0)
ax.set_ylim(np.min(u_test_min), np.max(u_test_max))
elif modelo == "onda":
u_test = u_sol_onda(x_vals, 0)
ax.set_ylim(np.min(u_test)*1.1, np.max(u_test)*1.1)
ax.set_xlim(0, L)
ax.grid(True)
if modelo == "calor":
t_anim = np.linspace(0, 500, 300)
else:
t_anim = np.linspace(0, 10, 300)
def init():
line.set_data([], [])
return line,
def update(frame):
t_val = t_anim[frame]
if modelo == "calor":
u_val = u_sol_calor(x_vals, t_val)
else:
u_val = u_sol_onda(x_vals, t_val)
line.set_data(x_vals, u_val)
ax.set_title(f"t = {t_val:.2f} s")
return line,
anim = FuncAnimation(
fig,
update,
frames=len(t_anim),
init_func=init,
interval=50,
blit=False
)
from IPython.display import HTML
HTML(anim.to_jshtml())
Out[23]:
BLOQUE 8 - Superficie 3D u(x,t)¶
Calor¶
In [11]:
from mpl_toolkits.mplot3d import Axes3D
# Mallado de espacio y tiempo
x_vals = np.linspace(0, L, 200)
t_vals = np.linspace(0, 500, 1000)
X, T = np.meshgrid(x_vals, t_vals)
# Matriz U(x,t)
U = np.zeros_like(X)
for i in range(len(t_vals)):
U[i, :] = u_sol_calor(x_vals, t_vals[i])
# Figura 3D
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
# Rotar la vista 3D
ax.view_init(elev=30, azim=45)
surf = ax.plot_surface(X, T, U, cmap='inferno', edgecolor='none')
ax.set_title("Evolución térmica u(x,t) — Ecuación de Calor")
ax.set_xlabel("Posición x")
ax.set_ylabel("Tiempo t")
ax.set_zlabel("Temperatura u")
fig.colorbar(surf, shrink=0.5, aspect=8)
plt.show()
Onda¶
In [24]:
from mpl_toolkits.mplot3d import Axes3D
# Mallado de espacio y tiempo
x_vals = np.linspace(0, L, 200)
t_vals = np.linspace(0, 80, 200) # la onda oscila rápido → tiempos pequeños
X, T = np.meshgrid(x_vals, t_vals)
# Matriz U(x,t)
U = np.zeros_like(X)
for i in range(len(t_vals)):
U[i, :] = u_sol_onda(x_vals, t_vals[i])
# Figura 3D
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
# Rotar la vista 3D
ax.view_init(elev=30, azim=45)
surf = ax.plot_surface(X, T, U, cmap='viridis', edgecolor='none')
ax.set_title("Evolución u(x,t) — Ecuación de Onda")
ax.set_xlabel("Posición x")
ax.set_ylabel("Tiempo t")
ax.set_zlabel("Desplazamiento u")
fig.colorbar(surf, shrink=0.5, aspect=8)
plt.show()
Bloque 9 - Solución con FTCS y comparaciones¶
Calor¶
In [14]:
def explicit_fd(f_init, L, a, T0, TL, nx, t_max):
dx = L/(nx-1)
# dt controlado por estabilidad CFL
dt = 0.45 * dx**2 / a
nt = int(t_max / dt) + 1
x_grid = np.linspace(0, L, nx)
times = np.linspace(0, t_max, nt)
u = f_init(x_grid)
u[0], u[-1] = T0, TL
hist = np.zeros((nt, nx))
hist[0] = u.copy()
r = a * dt / dx**2
for k in range(1, nt):
un = u.copy()
u[1:-1] = un[1:-1] + r*(un[2:] - 2*un[1:-1] + un[:-2])
u[0], u[-1] = T0, TL
hist[k] = u.copy()
return x_grid, times, hist
nx_num = 10
t_max = 80.0
print(f"Ejecutando FTCS con nx = {nx_num}")
f_init_fn = lambda X: u_sol_calor(X, 0)
x_num, t_num, U_num_hist = explicit_fd(
f_init_fn, L, alpha, T0, TL,
nx=nx_num,
t_max=t_max
)
U_analytic_hist = np.zeros_like(U_num_hist)
for i, tt in enumerate(t_num):
U_analytic_hist[i] = u_sol_calor(x_num, tt)
t_compare = [0, 1, 10, 40]
plt.figure(figsize=(9,6))
for tt in t_compare:
idx = np.argmin(np.abs(t_num - tt))
plt.plot(x_num, U_analytic_hist[idx], lw=2, label=f"AnalÃtica t={t_num[idx]:.2f}")
plt.plot(x_num, U_num_hist[idx], '--', lw=1.5, label=f"FTCS t={t_num[idx]:.2f}")
plt.title(f"Comparación AnalÃtica vs FTCS — nx={nx_num}")
plt.xlabel("x"); plt.ylabel("u(x,t)")
plt.grid(True); plt.legend()
plt.show()
x_surf, t_surf = x_num, t_num
X_surf, T_surf = np.meshgrid(x_surf, t_surf)
U_analytic_surf = U_analytic_hist
U_num_surf = U_num_hist
fig = plt.figure(figsize=(18,5))
# analÃtica
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X_surf, T_surf, U_analytic_surf, cmap='inferno')
ax1.set_title("AnalÃtica u(x,t)")
ax1.view_init(30,45)
# numérica
ax2 = fig.add_subplot(132, projection='3d')
ax2.plot_surface(X_surf, T_surf, U_num_surf, cmap='inferno')
ax2.set_title("FTCS u_num(x,t)")
ax2.view_init(30,45)
# error
ax3 = fig.add_subplot(133, projection='3d')
ax3.plot_surface(X_surf, T_surf, U_num_surf - U_analytic_surf, cmap='seismic')
ax3.set_title("Error u_num - u_analytic")
ax3.view_init(30,45)
plt.tight_layout()
plt.show()
Ejecutando FTCS con nx = 10
Onda¶
In [26]:
nx_num = 10
t_max = 10.0
f_init_fn = lambda X: u_sol_onda(X, 0)
x_num, t_num, U_num_hist = explicit_fd(
f_init_fn, L, a, T0, TL,
nx=nx_num,
t_max=t_max
)
U_analytic_hist = np.zeros_like(U_num_hist)
for i, tt in enumerate(t_num):
U_analytic_hist[i, :] = u_sol_onda(x_num, tt)
t_compare = [0.0, 1.0, 2.0, 4.0]
plt.figure(figsize=(9,6))
for tt in t_compare:
idx = np.argmin(np.abs(t_num - tt))
plt.plot(x_num, U_analytic_hist[idx], '-', lw=2, label=f'AnalÃtica t={t_num[idx]:.2f}')
plt.plot(x_num, U_num_hist[idx], '--', lw=1.5, label=f'FTCS t={t_num[idx]:.2f}')
plt.xlabel('x'); plt.ylabel('u(x,t)')
plt.title(f'Comparación AnalÃtica vs FTCS (onda) — nx={nx_num}')
plt.legend(); plt.grid(True)
plt.show()
x_surf, t_surf, U_num_surf = explicit_fd(f_init_fn, L, a, T0, TL, nx=nx_num, t_max=t_max)
U_analytic_surf = np.zeros_like(U_num_surf)
for i, tt in enumerate(t_surf):
U_analytic_surf[i, :] = u_sol_onda(x_surf, tt)
X_surf, T_surf = np.meshgrid(x_surf, t_surf)
fig = plt.figure(figsize=(18,5))
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X_surf, T_surf, U_analytic_surf, cmap='viridis', edgecolor='none')
ax1.set_title("AnalÃtica onda u(x,t)")
ax1.set_xlabel("x"); ax1.set_ylabel("t"); ax1.set_zlabel("u")
ax1.view_init(elev=30, azim=45)
ax2 = fig.add_subplot(132, projection='3d')
ax2.plot_surface(X_surf, T_surf, U_num_surf, cmap='viridis', edgecolor='none')
ax2.set_title("Numérica (FTCS) — NO exacta")
ax2.set_xlabel("x"); ax2.set_ylabel("t"); ax2.set_zlabel("u")
ax2.view_init(elev=30, azim=45)
ax3 = fig.add_subplot(133, projection='3d')
ax3.plot_surface(X_surf, T_surf, U_num_surf - U_analytic_surf, cmap='seismic', edgecolor='none')
ax3.set_title("Error FTCS - AnalÃtica")
ax3.set_xlabel("x"); ax3.set_ylabel("t"); ax3.set_zlabel("error")
ax3.view_init(elev=30, azim=45)
plt.tight_layout()
plt.show()
Bloque 10 - Análisis de error y convergencia¶
In [15]:
Error = U_num_hist - U_analytic_hist
L2 = np.sqrt(np.sum(Error**2, axis=1) / Error.shape[1])
Linf = np.max(np.abs(Error), axis=1)
plt.figure(figsize=(8,4))
plt.plot(t_num, L2, label='L2 (RMS) del error')
plt.plot(t_num, Linf, label='L∞ del error')
plt.xlabel('Tiempo t'); plt.ylabel('Error')
plt.title('Normas del error en el tiempo')
plt.legend(); plt.grid(True)
plt.show()
plt.figure(figsize=(10,5))
plt.contourf(x_num, t_num, Error, levels=60, cmap='seismic')
plt.colorbar(label='u_num - u_analytic')
plt.xlabel('x'); plt.ylabel('t')
plt.title('Mapa de calor del error E(x,t)')
plt.show()
t_final = min(40.0, t_num[-1])
nx_list = [50, 100, 200, 400]
errors_final = []
for nx_test in nx_list:
x_test, t_test, U_test = explicit_fd(f_init_fn, L, alpha, T0, TL, nx=nx_test, t_max=t_final)
U_num_final = U_test[-1, :]
U_analytic_final = u_sol_calor(x_test, t_final)
err_L2 = np.sqrt(np.sum((U_num_final - U_analytic_final)**2) / U_num_final.size)
errors_final.append(err_L2)
print("nx \t L2 error @ t_final")
for nxv, ev in zip(nx_list, errors_final):
print(f"{nxv}\t {ev:.6e}")
nx L2 error @ t_final 50 4.219587e-03 100 1.229874e-03 200 1.882337e-04 400 6.767316e-05
